function [ output ] = HousePriceElasticity( data , Params , options  )

st         = data.st ;
f          = data.f ;
W          = data.W ;
pop        = data.pop ;
u          = data.u ;
area       = data.area ;

houseprice = data.houseprice ; 

% recover zeta, v and bpv
run InvertZetaV.m

% House price contributions
I1 =   Params.omega .* u .* Params.b  ...
    .* W ./ bpv ./ ( 1 - Params.beta + Params.beta .* Ez ) ;
I2 = Params.omega .* ( 1 - u ) .* W ;
I3 = Params.psi .* W .* ( Ez ./ ( 1 - Params.beta + Params.beta .* Ez ) ) ;

% House prices

if options.HousePriceElasticity == 0
    eta = 0 ;
elseif options.HousePriceElasticity == 1
    
    % sample selection
    keepit = houseprice > 0 ;
    
    % LHS and RHS variables
    rhs = pop .* ( I1 + I2 + I3 ) ./ area ;
    RHS = log( rhs(keepit) ) ;
    LHS = log( houseprice(keepit) ) ;
    
    % uniform weights
    weights = ones(size(pop,1),1) ;
    weights = weights(keepit) ;
    
    % run regression
    OLS = fitlm( RHS , LHS , 'Weights' , weights ) ;
    
    % recover eta
    eta = 1 / OLS.Coefficients.Estimate(2) - 1 ;
    
end    
    
% Assign output
output.eta = eta ;


end

